#include <sys/types.h>
#include <sys/stat.h>

#include "Visualizator/visualizator.03.h"
#include "Vector/vector.03.h"
#include "Color/color.01.h"
#include "Timer/timer.01.h"
#include <vector>
#include <deque> 


const double kB = 1.3806488e-23;

using namespace std;

struct P {
	Vector<double> R;
	Vector<double> V;
	Vector<double> F;
};


struct PT {

	vector<P> p; 

	double m;
	double T;
	
	vector<double> e;
	vector<double> s; 
	vector<double> rc;
	
	Color color;
};



class MyVisualizator: public Visualizator {
	public:
		
		char **argv;
		
		vector<P> p;
		
		deque<PT> pt;
		
		bool coolingF;
		
		double rnd() {
			return (double)rand()/RAND_MAX;
		};
		
		
		void packToHCP(vector< P > &p, Vector<double> R, Vector<int> S, Vector<double> V, double t) {
			
			double dr = 1*0.977;
			
			grid_step = dr;
			
			for (int i=0; i<S.x; i++)
				for (int j=0; j<S.y; j++) 
					for (int k=0; k<S.z; k++){
					
						P pi;	
					
						pi.R.x = dr*i + (j%2==0 ? 0: 0.5*dr) + (k%2==0 ? 0: 0.5*dr);
						pi.R.y = 0.5*sqrt(3)*dr*j + (k%2==0 ? 0: sqrt(3.0)/6.0*dr);
						pi.R.z = sqrt(6.0)/3.0*dr*k;
					
						pi.V = Vector<double>(t*(1-2*rnd()), t*(1-2*rnd()), t*(1-2*rnd()));
					
						
						pi.R.x-= 0.5*S.x*dr;
						pi.R.y-= 0.5*S.y*0.5*sqrt(3)*dr;
						pi.R.z-= 0.5*S.z*sqrt(6.0)/3.0*dr;
						pi.R  += R;
						
						pi.V  += V;

						p.push_back(pi);	
					
				}
		}
		
		void packToHCP(PT &pt, Vector<double> R, Vector<int> S, Vector<double> V, double t) {
			
			double dr = pt.s[0]*0.9762;
			
			grid_step = dr;
			
			for (int i=0; i<S.x; i++)
				for (int j=0; j<S.y; j++) 
					for (int k=0; k<S.z; k++) {
					
						P pi;	
					
						pi.R.x = dr*i + (j%2==0 ? 0: 0.5*dr) + (k%2==0 ? 0: 0.5*dr);
						pi.R.y = 0.5*sqrt(3)*dr*j + (k%2==0 ? 0: sqrt(3.0)/6.0*dr);
						pi.R.z = sqrt(6.0)/3.0*dr*k;
					
						pi.V = Vector<double>(t*(1-2*rnd()), t*(1-2*rnd()), t*(1-2*rnd()));
						
						pi.R.x-= 0.5*S.x*dr;
						pi.R.y-= 0.5*S.y*0.5*sqrt(3)*dr;
						pi.R.z-= 0.5*S.z*sqrt(6.0)/3.0*dr;
						pi.R  += R;
						
						pi.V  += V;

						pt.p.push_back(pi);	
					
					}
		}
		
		MyVisualizator(){
			eye.R.z = eye.RM.z = 30;
			
			PT ptU;
			ptU.s.push_back(1);
			ptU.e.push_back(1);
			ptU.m = 1;
			packToHCP(ptU, Vector<double>(0,0,0), Vector<int>(32,32,32),  Vector<double>(0,0,0), 0);
			ptU.color = Color(1,0.2,0.1,1);
			pt.push_back(ptU);
			
			/*
			PT pt16;
			pt16.s.push_back(1);
			pt16.e.push_back(1);
			pt16.s.push_back(1);
			pt16.e.push_back(500);

			pt16.m = 1;
			packToHCP(pt16, Vector<double>(-10,0.5,0), Vector<int>(10,10,10),  Vector<double>(5,0,0), 0);
			pt16.color = Color(0.1,0.2,1,1);
			pt.push_back(pt16);
			*/
			
			perspectiveF = false;
			scale = -1;
			pointSize = 5;
			
			integrator.parent = this;
			integrator.start();
			coolingF = false;
			
			
			grid_step = 1;
		};
		
		class Integrator : public QThread {
			public:
				
				MyVisualizator *parent;

				typedef long unsigned int luint;
				
				dVector w;
				bool exitFlag;
				bool ssF, eF, pF, firstTimeF;
				double dt;
				double pot;
				double lU, eU, mU, tU;
				double etot, ek, ep, etot0, eerr;
				double temp;
				uint ntot;
				unsigned int step, nsteps;
				Timer tr, trf;
				double drij2max;
				unsigned int nstep;
				long unsigned int flops1, flops2, flops;
				
				Integrator():
				ssF(false),
				eF(false),
				pF(false),
				firstTimeF(true),
				nsteps(1<<32-1),
				flops1(0),flops2(0),flops(0)
				{
					lU = 1e-9; // m [1nm]
					eU = 1.60217646e-19; // J [1ev]
					mU = 1.660538921e-27; // kg [amu]
					tU = lU*sqrt(mU/eU); // sec
				}
				

				double find_lj_forces_(deque<PT> &pt) {
					
					luint loops(0);
					double pot = 0.0, f = 0.0;
					
					drij2max = SQR(3.0); 
					
					for (int i1=0;i1<pt.size();i1++) 
						for (int i2=0;i2<pt[i1].p.size();i2++) 
							memset(&pt[i1].p[i2].F, 0, sizeof(pt[i1].p[i2].F));
	
					
					for (int i1=0;i1<pt.size();i1++) 
					for (int i2=0;i2<pt[i1].p.size();i2++)
						for (int j1=0;j1<pt.size();j1++)
						for (int j2=0;j2<pt[j1].p.size();j2++) 
							if ( i1!=j1 || i2!=j2 ) { // 3+5 + 4+3+4+6+1 = 8 + 18 = 26
							
								
								double dx = pt[i1].p[i2].R.x - pt[j1].p[j2].R.x; // 1
								double dy = pt[i1].p[i2].R.y - pt[j1].p[j2].R.y; // 1
								double dz = pt[i1].p[i2].R.z - pt[j1].p[j2].R.z; // 1
								
								double r2 = dx*dx + dy*dy + dz*dz; // 5

								if (r2>drij2max) continue; 
							
								double r4 	= r2*r2; // 1
								double r6i 	= 1.0/(r2*r4); // 2
								double r12i = r6i*r6i;  // 1
								
								pot += r12i - 2*r6i;  // 3
								f = 12*r6i*r4*(r12i-r6i); // 4

								
								pt[i1].p[i2].F.x += dx*f; // 2
								pt[i1].p[i2].F.y += dy*f; // 2
								pt[i1].p[i2].F.z += dz*f; // 2
								
								loops++; //1 
							
							}

					
					
					luint n(0);
					for (int i=0;i<pt.size();i++) 
					n+=pt[i].p.size();
						
					flops1 = 8*(n*(n-1));
					flops2 = 18*loops;
					flops = flops1 + flops2;
					
					return pot;

				}					
				
				
				
				double find_lj_forces(deque<PT> &pt) {
					//printf("sfsf\n");
					luint loops(0);
					double pot = 0.0, f = 0.0;
					
					drij2max = SQR(3.0); 
					
					for (int i1=0;i1<pt.size();i1++) 
						for (int i2=0;i2<pt[i1].p.size();i2++) 
							memset(&pt[i1].p[i2].F, 0, sizeof(pt[i1].p[i2].F));
						
					 for (int i1=0;i1<pt.size();i1++) 
						for (int j1=0;j1<i1+1;j1++) {
							
							double eij = pt[i1].e[j1];
							
							for (int i2=0;i2<(j1==i1 ? pt[i1].p.size()-1:pt[i1].p.size());i2++) {
								P &pti1pi2 = pt[i1].p[i2];
								P *ptj1p = &pt[j1].p[0];
								for (int j2=(j1==i1 ? i2+1:0);j2<pt[j1].p.size();j2++) {
					
					
									double dx = pti1pi2.R.x - ptj1p[j2].R.x; // 1
									double dy = pti1pi2.R.y - ptj1p[j2].R.y; // 1
									double dz = pti1pi2.R.z - ptj1p[j2].R.z; // 1
								
									double r2 = dx*dx + dy*dy + dz*dz; // 5
						
									if (r2>drij2max) continue; 
							
									double r4 	= r2*r2; // 1
									double r6i 	= 1.0/(r2*r4); // 2
									double r12i = r6i*r6i;  // 1
								
									pot += eij*(r12i - 2*r6i);  // 3
									f = eij*12*r6i*r4*(r12i-r6i); // 4
								
									double dxf = dx*f; // 1
									double dyf = dy*f; // 1
									double dzf = dz*f; // 1
									
									pti1pi2.F.x += dxf; // 1
									pti1pi2.F.y += dyf; // 1
									pti1pi2.F.z += dzf; // 1
								
									ptj1p[j2].F.x -= dxf; // 1
									ptj1p[j2].F.y -= dyf; // 1
									ptj1p[j2].F.z -= dzf; // 1
					
									loops++; //1 
								}
							}
						}
				
					luint n(0);
					for (int i=0;i<pt.size();i++) 
					n+=pt[i].p.size();
				
					flops1 = 8*(n*(n-1)/2);
					flops2 = 21*loops;
					flops = flops1 + flops2;
					
					return pot;
				}
				
				
				
				double ___find_lj_forces(deque<PT> &pt) {
					
					luint loops(0);
					double pot = 0.0, f = 0.0;
					
					drij2max = SQR(3.0); 
					
					for (int i1=0;i1<pt.size();i1++) 
						for (int i2=0;i2<pt[i1].p.size();i2++) 
							memset(&pt[i1].p[i2].F, 0, sizeof(pt[i1].p[i2].F));
	

					for (int i1=0;i1<pt.size();i1++) 
					for (int i2=0;i2<pt[i1].p.size();i2++) {
						P &pi = pt[i1].p[i2];
						
						for (int j1=0;j1<pt.size();j1++) {
						P *pj1 = &pt[j1].p[0];
						
						for (int j2=0;j2<pt[j1].p.size();j2++) 
							if ( i1!=j1 || i2!=j2 ) { // 3+5 + 4+3+4+6+1 = 8 + 18 = 26

								double dx = pi.R.x - pj1[j2].R.x; // 1
								double dy = pi.R.y - pj1[j2].R.y; // 1
								double dz = pi.R.z - pj1[j2].R.z; // 1
								
								double r2 = dx*dx + dy*dy + dz*dz; // 5
						
								if (r2>drij2max) continue; 
							
								double r4 	= r2*r2; // 1
								double r6i 	= 1.0/(r2*r4); // 2
								double r12i = r6i*r6i;  // 1
								
								pot += (r12i - 2*r6i);  // 3
								f = 12*r6i*r4*(r12i-r6i); // 4
								
								
								
								pi.F.x += dx*f; // 2
								pi.F.y += dy*f; // 2
								pi.F.z += dz*f; // 2
								
								loops++; //1 
							
							}
						}

					}
					
					luint n(0);
					for (int i=0;i<pt.size();i++) 
					n+=pt[i].p.size();
						
					flops1 = 8*(n*(n-1));
					flops2 = 18*loops;
					
					
					flops = flops1 + flops2;
					
					return pot;

				}	

				
				double find_lj_forces(vector<P> &p) {
					
					luint loops(0);
					double pot = 0.0, f = 0.0;
					
					drij2max = SQR(3.0); 
					
					for (int i=0;i<p.size();i++) memset(&p[i].F, 0, sizeof(p[i].F));
	
						
					for (int i=0;i<p.size()-1;i++) 
					for (int j=i+1;j<p.size();j++) { // 3 + 5  +  1+2+1 + 2+6 + 1+1+1 + 1+1+1+1+1+1 = 8 + 21 =  29
						
						double dx = p[i].R.x - p[j].R.x; // 1
						double dy = p[i].R.y - p[j].R.y; // 1
						double dz = p[i].R.z - p[j].R.z; // 1
        
						double r2 = dx*dx + dy*dy + dz*dz; // 5
			
						if (r2>drij2max) continue; 
         
						double r4 	= r2*r2; // 1
						double r6i 	= 1.0/(r2*r4); // 2
						double r12i = r6i*r6i;  // 1
        
						pot += r12i - 2*r6i;  // 2
						f = 12*r6i*r4*(r12i-r6i); // 4
		  
						double dxf = dx*f; // 1
						double dyf = dy*f; // 1
						double dzf = dz*f; // 1
		  	
						p[i].F.x += dxf; // 1
						p[j].F.x -= dxf; // 1
						p[i].F.y += dyf; // 1
						p[j].F.y -= dyf; // 1
						p[i].F.z += dzf; // 1
						p[j].F.z -= dzf; // 1
						
						loops++; //1 
						
					}
					
					flops1 = 8*((luint)p.size()*(p.size()-1))/2;
					flops2 = 21*loops;
					flops = flops1 + flops2;
					return pot;
				}	
				
				
				
				void run(){
					
					
					exitFlag = false;
					
					deque<PT> &pt = parent->pt;
					
					step=0;
					
					etot0 = 0;
					for(int i1=0; i1<pt.size(); i1++)
					for(int i2=0; i2<pt[i1].p.size(); i2++)
					etot0+= pt[i1].p[i2].V.sqr(); 
					
					etot0*=0.5;
					etot0+= find_lj_forces(pt);
					
					
					
					ntot=0;
					for(int i=0; i<pt.size(); i++) ntot+=pt[i].p.size();
					
					dt = 1e-10*pow(ntot, 2); if (dt>0.01) dt = 0.01;
					
					do {

						tr.set();
						
						for(int i1=0; i1<pt.size(); i1++)
						for(int i2=0; i2<pt[i1].p.size(); i2++) {
							pt[i1].p[i2].R+= (pt[i1].p[i2].V + .5*pt[i1].p[i2].F*dt)*dt;
							pt[i1].p[i2].V+= .5*pt[i1].p[i2].F*dt;
						}
						
						trf.set();
 						ep = find_lj_forces(pt);
						trf.get();
						
						for(int i1=0; i1<pt.size(); i1++)
						for(int i2=0; i2<pt[i1].p.size(); i2++)
							pt[i1].p[i2].V+= .5*pt[i1].p[i2].F*dt;
						/*
						if (parent->coolingF) 
						for(int i1=0; i1<pt.size(); i1++)
						for(int i2=0; i2<pt[i1].p.size(); i2++)
							pt[i1].p[i2].V*= 0.999; 
						*/
						
						ek=0;
						for(int i1=0; i1<pt.size(); i1++)
						for(int i2=0; i2<pt[i1].p.size(); i2++)
						ek+= pt[i1].p[i2].V.sqr(); ek*=.5;
						
						
						for(int i1=0; i1<pt.size(); i1++) {
							
							Vector<double> pti1Vcm;
							double pti1ek(0),pti1eki(0),pti1ekc(0); 
							
							for (int i2=0; i2<pt[i1].p.size(); i2++) {
								pti1ek+=pt[i1].p[i2].V.sqr();
								pti1Vcm+=pt[i1].p[i2].V; 
							}
							pti1ek*=.5;
							pti1ekc=.5*pti1Vcm.sqr()/pt[i1].p.size();
							pti1eki=pti1ek-pti1ekc;
							pt[i1].T = (2.0/3)*pti1eki/pt[i1].p.size();
						}
						
						if (step%20==0) printf("% 5i % .4f  % .4e [s]    % .4f % .4f [K]\n", step, step*dt, step*dt*tU, pt[0].T, pt[0].T*eU/kB);
						
						etot=ek+ep;	
						eerr= (etot-etot0)/etot0;
						
						step++;

						tr.get();
						
					} while (step<nsteps && !exitFlag);
					
				}
				
				
				
				
				
				void _run(){
					
					
					exitFlag = false;
					
					vector<P> &p = parent->p;	
					
					step=0;
					
					etot0 = 0;
					for(int i=0; i<p.size(); i++)
					etot0+= p[i].V.sqr(); etot0*=0.5;
					etot0+= find_lj_forces(p);
					
					dt = 1e-9*pow(p.size(), 2); if (dt>0.01) dt = 0.01;					

					do {

						tr.set();
						
						for(int i=0; i<p.size(); i++) {
							p[i].R+= (p[i].V + .5*p[i].F*dt)*dt;
							p[i].V+= .5*p[i].F*dt;
						}
						trf.set();
 						ep = find_lj_forces(p);
						trf.get();
						
						for(int i=0; i<p.size(); i++)
							p[i].V+= .5*p[i].F*dt;
							
						for(int i=0; i<p.size(); i++)
						if (parent->coolingF) p[i].V*= 0.999; 
						
						ek=0;
						for(int i=0; i<p.size(); i++)
						ek+= p[i].V.sqr(); ek*=.5;
						
						
						
						
						
						etot=ek+ep;	
						eerr= (etot-etot0)/etot0;
						
						step++;
						
						/*
						
						if (ssF && step%5==0) {
							
							char filename[256];
							sprintf(filename, "ss/%03i.txt",step);
							
							FILE *out = fopen(filename, "w");
							for(int i=0; i<p.size(); i++)
							fprintf(out, "% .4e % .4e % .4e\n", p[i].R.x, p[i].R.y, p[i].R.z);
							fclose(out);
							
						
						}
						
						
						if (eF && step%1==0) {
							
							char filename[256];
							sprintf(filename, "%s.e.txt", parent->argv[0]);
							
							FILE *out = fopen(filename, (firstTimeF ? "w" : "a"));
							fprintf(out,"%6i %.7e  %.7e %.7e %.7e %.7e\n", step, step*dt, ek, ep, e, e_err);
							fclose(out);
							
							firstTimeF = false;	
						}
						
						
						if (pF && step%10==0) {
							
							char filename[256];
							sprintf(filename, "%s.p.txt", parent->argv[0]);
							
							FILE *out = fopen(filename, (firstTimeF ? "w" : "a"));
							fprintf(out,"%6i %15lu %15lu  %15lu  %8.3f %8.4f\n", step, flops1,flops2,flops, flops/trf.sin*1e-9, 1/tr.sin);
							fclose(out);
							
							firstTimeF = false;	
						}
						
						*/
						tr.get();
						
						
					} while (step<nsteps && !exitFlag);
					
				}
				
				~Integrator(){ exitFlag = true; }
				
		} integrator;
		
		
		virtual void printInfo(){
			
			Visualizator::printInfo();
			
			
			char str[512];
			sprintf(str, 
			"ek: %.2e\nep: %+.2e\ne: %+.2e\ne_err: %+.2e\nT:%.2f\n\nstep: % .2e\ndt: %.2e\nt: %.2f\nsteps/s: %g\nsteps/s: %g (force)\n\np.size(): %i\nflops1: %lu\nflops2: %lu\nflops:   %lu\nratio: %.2f\nGlops/s: %.3f\n", 
					  integrator.ek, integrator.ep, integrator.etot, integrator.eerr, pt[0].T,(double)integrator.step, integrator.dt, integrator.step*integrator.dt,
					  1.0/integrator.tr.sin, 1.0/integrator.trf.sin, integrator.ntot,
					  integrator.flops1, integrator.flops2, integrator.flops, (double)integrator.flops1/integrator.flops2, (double)integrator.flops/integrator.trf.sin*1e-9
					  );
			
			glColor3f (.4,.4,.4);
			glPrint(10,85, str, GLUT_BITMAP_HELVETICA_12);
			
		}
		
		void paintGL() {
			
			Visualizator::paintGL();
			
			
			for(int k=0; k<pt.size(); k++) {
				glBindBuffer(GL_ARRAY_BUFFER_ARB, k+1);
				glBufferData(GL_ARRAY_BUFFER_ARB, pt[k].p.size()*sizeof(pt[k].p[0]), &(pt[k].p[0].R), GL_STATIC_DRAW_ARB);
				drawPointArray(k+1, pt[k].p.size(), pt[k].color, pt[k].color, pointSize, GL_ONE, GL_ONE, false, program, texture, sizeof(p[0]), 0, GL_DOUBLE);
			}
			
		};
		
		void keyPressEvent(QKeyEvent *event){
			Visualizator::keyPressEvent(event);
			
			if (event->key() == Qt::Key_PageUp)	scale+= -1.0/scale_steps;
			if (event->key() == Qt::Key_PageDown) scale+= +1.0/scale_steps;
			
			if (event->key() == Qt::Key_C) coolingF=!coolingF;
			
			
		}
	
};



int main(int argc, char *argv[]) {
	
	glutInit(&argc, argv); //initialize glut
	QApplication app(argc, argv);
	
	MyVisualizator visualizator;
	
	bool windowF(true);
	visualizator.argv=argv;
	
	for (int i=1;i<argc;i++)
		if (!strcmp(argv[i],"-nw")) windowF=false;
		else if (!strcmp(argv[i],"-ssF")) { visualizator.integrator.ssF=true;  mkdir("ss", 0777);  }
		else if (!strcmp(argv[i],"-eF")) visualizator.integrator.eF=true;
		else if (!strcmp(argv[i],"-pF")) visualizator.integrator.pF=true;
		else if (!strcmp(argv[i],"-nsteps")) visualizator.integrator.nsteps = atoi(argv[++i]);
	
	
	visualizator.resize(800, 800);
	visualizator.setWindowTitle( QString("Visualizator Lite") );
	if (windowF) visualizator.show();
	
	return app.exec();
	
}
